#####################################################################################################
# This script generates a series of maps.                                                           #
# It draws on existing data series on the location of lakes, rivers, canals, and railroad,          #
# as well as the topography of the state. Only the first of these maps is included in the main      #
# text.                                                                                             #
#                                                                                                   #
# (1) County and town level proportion of enfranchised (assembly) adult men in NYS in 1814.         #
# (2) County and town level proportion of enfranchised (governor/senate) adult men in NYS in 1814.  # 
#                                                                                                   #
#####################################################################################################

### Start by loading packages used across different maps. 
rm(list = ls())

library(shapefiles)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(ggplot2)
library(readxl)
library(dismo)
library(elevatr)
library(raster)
library(rgeos)
library(ggmap) 
library(sf) 

### Provide version list. 
sessionInfo()

### Set the API key to get the Google geocoded location of cities. For users without a key, longitude and latitude are provided in one of the input files.
#register_google('[API key here]')
dir <-paste("../", sep="")
setwd(dir)

### Create a heat map function that will be used across different maps. 

## Create a heat map function
plot.heat <- function(states.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(states.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(states.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  states.map@data$zCat <- cut(states.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(states.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(states.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(states.map@data$zCat) <- cutpointsColors
  plot(states.map,border=gray(.8), lwd=bw,axes = FALSE, las = 1,col=as.character(states.map@data$zCat))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}

### Define color scheme and create index for figures
color.palette1 =  colorRampPalette(c("black","white"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(6)
index = c(cols[1], cols[3], cols[4], cols[5], cols[6])

################################################################################
###         New York - 1814 Assembly Voters and 1814 Governor Voters         ###
################################################################################

## Set year interval to make some of the shapefiles more manageable. 
year <- paste("1810")
year2<- paste("1820")

## To generate Assembly figure (main text) make office = "Assembly"
## To generate Governor figure (supplemental) make office = "Governor"
office = "Governor"
if(office=="Governor"){
	office1 	<-paste("govPercent")
	disf		<-paste("totDisfGov")
	label  	<-paste("Figures/Supplemental/Supplemental-Figure-1.pdf", sep="")
	pt.cex.size = c(6.5, 3, 1, .5,0.01)
	pt.adj	= 200
	lg1 		= "1300"
	lg2		= "600"
	lg3		= "200"
	lg4		= "100"
	lg5		= "2"
	lgb1		= "10-25%"
	lgb2		= "25-35%"
	lgb3		= "35-45%"
	lgb4		= "45-50%"
	lgb5		= "50-55%"
	brk		= c(10,25,35,45,50,55)
} else if(office=="Assembly"){
	office1 = paste("legPercent")
	disf		= paste("totDisf")
	label  	<-paste("Figures/Figure-2.pdf", sep="")
	pt.cex.size = c(6.5, 3, 1, .5,0.01)
	pt.adj	= 100
	lg1 		= "650"
	lg2		= "300"
	lg3		= "100"
	lg4		= "50"
	lg5		= "1"
	lgb1		= "30-60%"
	lgb2		= "60-65%"
	lgb3		= "65-70%"
	lgb4		= "70-80%"
	lgb5		= "80-100%"
	brk		= c(30,60,65,70,80,100)
}


## Create list of inputs for shape/voter files. Three shapefiles: NYS, NYC, Albany, and Schenectady.
## Two voter files, one at the county level (% enfranchised) and one at the town level (number of
## adult men disenfranchised). The voter files are based on the New York State Censuses. 

input1 <- paste("ShapeFiles/NYS1814.shp", sep="")
input2 <- paste("Data/nycensus1814.csv", sep="")
input3 <- paste("ShapeFiles/NYC.shp", sep="")
input4 <- paste("ShapeFiles/Albany.shp", sep="")
input5 <- paste("Data/nycounties1814.csv", sep="")
input6 <- paste("ShapeFiles/Schenectady.shp", sep="")

## Load rivers, canals, the NYS counties, and elevation data. Lakes and canals can be included
## but I have commented them out here, since the Erie Canal was not yet completed and the lakes
## can sometimes cover over the whole plot. 

rivers <- readOGR("Shapefiles/Lakes_and_Rivers_Shapefile_NA_Lakes_and_Rivers_data_hydrography_l_rivers_v2.shp", layer = "Lakes_and_Rivers_Shapefile_NA_Lakes_and_Rivers_data_hydrography_l_rivers_v2")
rivers <- rivers[rivers@data[["COUNTRY"]]=="USA" | rivers@data[["COUNTRY"]]=="CAN",]
rivers <- spTransform(rivers, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
#lakes <- readOGR("Shapefiles/Lakes_and_Rivers_Shapefile_NA_Lakes_and_Rivers_data_hydrography_p_lakes_v2.shp", layer = "Lakes_and_Rivers_Shapefile_NA_Lakes_and_Rivers_data_hydrography_p_lakes_v2")
#canals <- readOGR("Shapefiles/19thC_Canals_March2017", layer = "19thC_Canals_March2017")
#canals <- canals[which(canals$OPENED < year & canals$CLOSED > year),]
#canals <- spTransform(canals, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
counties.map <- spTransform(counties.map, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

### Elevation is included to get a sense of whether enfranchisement varied with topography.
### Earlier versions of "elevatr" package allowed for use of sp objects. Version 0.99.0 requires sf or terra.
### Quick workaround without rewriting code is to convert counties.map to sf and then reload it as sp. 
### get_elev_raster requires internet connection to access topographic info.
### As noted below when plotting figure, elevation is not needed for reproduction and so this section 
### can be commented out. 
counties.map <- st_as_sf(counties.map)
elevation    <- get_elev_raster(counties.map, z = 8, prj=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
counties.map <- spTransform(counties.map, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

## Load NYC/Albany/Schenectady shape files and voter data

NYC.map <- readShapeSpatial(input3, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
Albany.map <- readShapeSpatial(input4, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
Schenectady.map <- readShapeSpatial(input6, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

data <- read.csv(input2)
data.counties <- read.csv(input5)

data.clean <- data[complete.cases(data$latitude, data$longitude),]
data.clean1 <-data.clean[which(data.clean$city==0), ]
data.clean2 <-data.clean[which(data.clean$city==1), ]
myvars <- c(disf, "longitude", "latitude")
data.clean1 <-data.clean1[myvars]
data.clean2 <-data.clean2[myvars]
names(data.clean1 )[1]<-"totDisf"
names(data.clean2 )[1]<-"totDisf"

names(data.counties )[28]<-"FIPS"
## Need to be dropping Warren county
data.counties <- subset(data.counties, FIPS!="36133")

counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,data.counties, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)==office1]<-"noise"

NYC.map@data$order <- 1:nrow(NYC.map)
NYC.map@data <- merge(NYC.map@data,data.counties, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
NYC.map@data<-NYC.map@data[order(counties.map@data$order),]
names(NYC.map@data)[ names(NYC.map@data)==office1]<-"noise"

Albany.map@data$order <- 1:nrow(Albany.map)
data.albany <- data.counties
Albany.map@data <- merge(Albany.map@data,data.albany, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
Albany.map@data<-Albany.map@data[order(counties.map@data$order),]
names(Albany.map@data)[ names(Albany.map@data)==office1]<-"noise"

Schenectady.map@data$order <- 1:nrow(Schenectady.map)
data.albany <- data.counties
Schenectady.map@data <- merge(Schenectady.map@data,data.albany, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
Schenectady.map@data<-Schenectady.map@data[order(counties.map@data$order),]
names(Schenectady.map@data)[ names(Schenectady.map@data)==office1]<-"noise"

## Create Figure 2 or Supplemental-Figure-1
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(0,0,0,0), mar=c(0,0,0,0))

plot(counties.map,density=c(30))
par(new=TRUE)
plot.heat(counties.map,z="noise",breaks=brk, plot.legend=FALSE)
plot(rivers, col = "dodgerblue", lwd = 1, add = TRUE)
#plot(canals, col = "blue", lwd = 1, add = TRUE)

### Create the topographical shading. This is not essential to the replication. 
### If user has difficulty getting elevation data from get_elev_raster above, then
### section can be commented out.
mosaic <- elevation *10
slope <- terrain(mosaic, opt="slope", unit='radians')
aspect <- terrain(mosaic, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)
hillshade2 <- aggregate(hillshade , fact = 5 , method = "bilinear" )
hillshade2 <- focal(hillshade2, w=matrix(1/9, nc=3, nr=3), mean)
slope2 <- aggregate(slope , fact = 5 , method = "bilinear" )
slope2 <- focal(slope2, w=matrix(1/9, nc=3, nr=3), mean)
plot(hillshade2, col=grey.colors(100, start=0, end=1), legend=F, alpha=.2, add=TRUE)
plot(slope2, col=grey.colors(100, start=1, end=0), legend=F, alpha=.2, add=TRUE)

### Plot location of towns / disenfranchised adult free men.
points(data.clean1$longitude, data.clean1$latitude, cex = (as.numeric(data.clean1$totDisf)+1)/pt.adj, pch=21,  bg=alpha(rgb(1,1,1), 0.6),col="black", lwd=.4)

### Insert the two indices, for counties and towns.
legend('bottomleft', fill=index , ncol=1,
       x.intersp = 1, xjust=0, yjust=0,cex=.9, bg="white", box.col = "white",
       legend = c(paste(lgb1, "Adult Free Men\nEnfranchised"), lgb2, lgb3, lgb4, lgb5))

legend(x = c(-78.1, -76.1), y = c(40, 41.4), y.intersp= 1.07, pch=21, cex=.9, 
       legend=c(paste(lg1, "Adult Free Men\nDisenfranchised"), lg2, lg3, lg4),
       pt.cex=pt.cex.size, x.intersp = 2, bg="white", box.col = "white")

### Insert the map insets for NYC/Albany/Schenectday
par(fig = c(0.02, .17, 0.65, .9),bg="white", new = T)
plot.heat(NYC.map,z="noise",breaks=brk, plot.legend=FALSE)
points(data.clean2$longitude, data.clean2$latitude, cex = (as.numeric(data.clean2$totDisf)+1)/(pt.adj*2), pch=21,  bg=alpha(rgb(1,1,1), 0.6),col="black", lwd=.4)
mtext("New York City\n and Brooklyn", cex=.7, adj=.5, family="sans")
box(lty = 'solid', col = 'black')

par(fig = c(.28, .4, 0.7, .9),bg="white", new = T)
plot.heat(Albany.map,z="noise",breaks=brk, plot.legend=FALSE)
points(data.clean2$longitude, data.clean2$latitude, cex = (as.numeric(data.clean2$totDisf)+1)/pt.adj, pch=21,  bg=alpha(rgb(1,1,1), 0.6),col="black", lwd=.4)
mtext("Albany", cex=.7, adj=.5, family="sans")
box(lty = 'solid', col = 'black')

par(fig = c(0.85, .95, 0.75, .95), bg="white", new = T)
plot.heat(Schenectady.map,z="noise",breaks=brk, plot.legend=FALSE)
points(data.clean2$longitude, data.clean2$latitude, cex = (as.numeric(data.clean2$totDisf)+1)/pt.adj, pch=21,  bg=alpha(rgb(1,1,1), 0.6),col="black", lwd=.4)
mtext("Schenectady", cex=.7, adj=.5, family="sans")
box(lty = 'solid', col = 'black')

dev.off()
